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ABSTRACT 

This thesis analyzes the problem of motion stability of submarines in free positive 
buoyancy ascent under casualty conditions such as control surface jam and loss of 
propulsion system response. We employ fully nonlinear, coupled six degree of freedom 
equations of motion and we allow response to occur in combined vertical and horizontal 
planes. Continuation and homotopy theory techniques are utilized to trace all possible 
steady state solutions in six degrees of freedom, while local perturbation reveals their 
stability properties. Vehicle geometric properties and control surface deflections are used 
as primary bifurcation parameters. Regions in parameter spaces are identified where 
extreme sensitivity of solutions to geometric properties and hydrodynamic modeling is 


present. 
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I. INTRODUCTION 


The dynamic response of a submarine under casualty conditions constitutes a 
crucial, and frequently limiting, factor in establishing the vehicle’s submerged operating 
envelope. As basic casualty conditions in the context of this work, we mean the loss of 
control surface and/or propulsion system response, and a flooding casualty where the 
boat must either be brought .to the surface or stabilized to a new operating depth. Of 
particular significance is the study of the ability of the boat to recover from a control 
surface jam. There exist several factors which determine the severity of such a situation 
as well as the recovery procedures; the initial conditions (speed, depth, pitch angle, etc.) 
during the jam, the actual control surface angles, reversing time and backing power of 
the propulsion system, the ability to blow ballast, and the time between recognition of 
casualty and initiation of proper recovery procedures. To this end, it is crucial that we 
have a clear understanding of the dynamics of the boat during an emergency situation. 

The traditional methods for establishing dynamic stability of motion concentrate 
mainly on eigenvalue analysis during small perturbations around nominal straight line 
paths Clayton and Bishop [Ref. 1]. Two indices are utilized, a stability index G, for the 
vertical plane and G, for the horizontal plane Roddy [Ref. 2]. In terms of the slow 


motion derivatives, these indices are given by: 


_ M,(Z,+m) 
” Zale) 


Nglt¥.—m) 
Yee Xam) 


Positive values for these indices, whose usage dates back to the 1950’s, indicates motion 
stability in the corresponding plane. The underlying assumption in these criteria is that 
during normal straight line motions, the coupling between horizontal and vertical plane 
motions is relatively weak, and can be neglected. Although vortex shedding and flow 
separation introduces a certain degree of coupling, the above assumption has been proven 
quite useful in design and analysis. However, for a high-speed fast-maneuvering 
submarine operating at the extremes of her submerged operating envelope, or during 
emergency situations, the above assumption of uncoupled motions breaks down. High 
amplitude motions may take place in all six degrees of freedom, and the nonlinear 
interactions between the various modes of motion become more pronounced. Therefore, 
we have to carefully consider the motion characteristics allowing for coupling between 
horizontal and vertical planes. 

The implications of nonlinear effects and coupling are numerous. In the case of roll 
motion, which is one of the most critical responses, there is growing evidence of 
complicated dynamics and chaotic response under certain excitations Falzarano, Shaw and 


Troesch [Ref. 3], Taz Ul Mulk and Falzarano [Ref. 4]. 


The motion of a submarine in the ocean environment involves some of the most 
complex fluid-structure interaction problems. For example, the steady and unsteady 
characteristics of the stern, and especially with regard to the effect of the wake and 
rolled-up body vortices on propulsor unsteadiness, and consequences for noise and 
hull vibrations are important problems of great concern. At present, too little is 
known in depth about the flow at the stern of ships and submarines and the 
footprints of their ensuing wakes in homogenous and stratified media Sarpkaya 
[Ref. 5]. 
Typical vortex structure about a submerged body can be seen in Figure 1, which is taken 


from Lugt [Ref. 6]. 


Ca ' 
—SSS ae 


Figure 1. Typical Vortex Structure About a Submerged Body (From Lugt, 1981). 


When an axisymmetric body moves at a sufficiently large angle of attack, the 
boundary layer vorticity may lift off the surface and sheets of vorticity are 


convected downstream while rolling up into streamwise vortices on the leeward 


side of the body. This separation and roll-up phenomenon, known as the crossplane 
separation, occurs on almost all maneuvering bodies at sufficiently high angles of 
attack. The roll-up of the vortices and their subsequent evolution are extremely 
important for the proper design of a submerged body and more important for the 
determination of the interaction of the vortices with the stern and the propulsion 
system. It is this interaction that determines the maneuverability and control of a 
submerged body Sarpkaya [Ref. 5]. 
In the case of submarine motions, there is evidence of bifurcation phenomena and 
extreme sensitivity of response to initial conditions and control actions during emergency 


TYPICAL SIMULATION 
20 . 
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Figure 2. Typical Simulation Results. 





ascent scenarios such as recovery from a dive plane jam. As motivation for the analysis 
that follows we present typical simulation results in Figure 2. The time simulation is in 
terms of vehicle roll angle ® versus time for 2% excess buoyancy with the buoyancy 
force located 1% of the vehicle length forward of the center of gravity, -6 degrees dive 
plane angle, and 0.1 feet metacentric height. Four different recovery actions are shown, 
all parametrized by the applied rudder angle in degrees. It can be seen that although 
zero rudder angle appears to bring the vehicle roll angle back to zero in the shortest 
time, the response does not persist. Small nonzero values for the rudder angle develop 
excessive roll angles in a divergent way, while larger rudder angles reduce the amount 
of roll. Clearly the vehicle effective rudder angle is largely affected by the amount of 
vorticity and currents in the flow field, and its actual value can not be known exactly. 
Since situations like the one presented in the graph should be avoided, we need to 
develop a mechanism for assessing those regions in the parameter space where the 
response can not be simulated with confidence. 

Stability in emergency ascent has been studied in the vertical plane. In Booth [Ref. 
7] and [Ref. 8], the vehicle response was distinguished into either a nearly vertical ascent 
or predominantly forward motion. In Papoulias and McKinley [Ref. 9], it was found that 
the above distinction is not always meaningful as a result of the many parameters that 
affect the problem. This latter study maintained some vertical plane restrictions, although 
it indicated the potential existence of pitchfork bifurcations which led to coupled out-of- 
plane solutions. In this work, we relax the requirement for vertical plane motions, and 


we analyze the stability properties of all possible steady states in six degrees of freedom. 


We employ a combination of singularity theory Golubitsky and Schaeffer [Ref. 10], 
bifurcation theory Guckenheimer and Holmes [Ref. 11], and numerical continuation 
methods Seydel [Ref. 12] in order to capture all steady state solutions that are physically 
admitted by the coupled nonlinear equations of motion. The primary bifurcation 
parameters used are: amount and location of excess buoyancy, diveplane and rudder 
deflection, and the metacentric height. Solution branching is shown to occur in various 
forms, including single and multiply connected pitchfork bifurcations, separation of 
solutions, hysteresis, and teardrop branches. Dynamic loss of stability in the form of 
Poincare-Andronov-Hopf bifurcations to periodic solutions is also identified. It is shown 
that for certain ranges of parameters, these periodic solutions persist in the vicinity of 
inverted pendulum steady states. Finally, we summarize our results in the form of 
bifurcation graphs which identify parameter regions with qualitatively different 
asymptotic response characteristics. For demonstration purposes, all computations in this 
work are performed for the Swimmer Delivery Vehicle, a 17.4 feet vehicle, for which 
a complete set of hydrodynamic and geometric properties is available Smith, Crane and 
Summey [Ref. 13]. Unless otherwise specified, all results in this work are presented in 
dimensional form, linear dimensions in feet, velocities in ft/sec, angular deflections in 


degrees, and time in seconds. 





Ii. PROBLEM FORMULATION 


A. PROBLEM STATEMENT 

Controlling emergency ascent situations on submersible vehicles such as dive plane 
jam recovery is of concern to the submariners. In order to control such situations, one 
must first be able to predict the dynamic response of positively buoyant submersibles. 

Dynamic response equations of motion describe the maneuvering characteristics of 
submersible vehicles for six degrees of freedom. These equations assume constant 
coefficients for hydrodynamic forces and moments approximated by zero frequency added 
mass and damping terms plus the quadratic terms for drag forces. The constant 
coefficients vary for each vehicle and depend on such things as vehicle body shape, 
location and magnitude of vehicle buoyancy, position of bow and stern planes, position 
of rudder, vehicle speed, vehicle mass characteristics, vehicle hydrodynamic coefficients, 
propeller rpm and control surface inputs. This thesis uses the equations of motion and 
hydrodynamic coefficients for a submerged Mark [IX Swimmer Delivery vehicle (SDV) 
developed by Smith, Crane and Summey [Ref. 13] to forecast the dynamic behavior of 
submersibles in a positively buoyant condition. 

For submarines where high amplitude motions may take place in all six degrees of 
freedom, nonlinear interactions between the various modes of motion may become more 
pronounced. In particular, there is growing evidence of bifurcation phenomena during 


high speed maneuvering and emergency ascent scenarios such as recovery from a dive 


plane jam. In these cases, it is no longer true that decoupled linear analysis techniques 
are sufficient and one is forced to consider the true character of six degrees of freedom 
motions. So, in this thesis, six degrees of freedom equations of motions have been used 
to investigate branching behaviors of the equations of motions. This is done under the 


condition of changing certain system parameters while keeping others constant. 


B. EQUATIONS OF MOTION 

The equations of motion given below are referenced to a nght-hand orthogonal axis 
system fixed in the body as shown in Figure 3. The origin of this body axis system is 
located in the vehicle’s X-Z plane of symmetry on the mid-body centerline behind the 
vehicle’s nose. Positive directions for control surface deflections, forces and angular 
moments, and linear and angular velocity components are shown in Figure 3. In the 
following equations, differentiation with respect to time is denoted by a dot over a 
quantity. The six degrees of freedom equations of motion for a submarine in surge, 


sway,heave, roll, pitch, and yaw, respectively, are: 


mx (u-vr+wq-xX,(q*+r7) +y,(pq-r) +2,.(pr+q) ] =X,+X,+X. ae 
inx [V+ur-wp+x,(pq+tr) -y, (p*+r*)4z,(qr=-p) )=Y,,+Y,+Y. ae 
mx [W-uq+vp+X,_(pr-q) +¥g(qr+p) -Zg(p?+q*)] =Zy+Zy+Z_ (3) 


I,p+ (I,-1,) quad, (pr=@) =f plee-= *) = aieg es) (4) 
+mx ly, (w-uqg+vp) -Z,(v+ur-wp) ] =K,+K,+Ko 
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yer (I,-I,) pr-I,,,(qr+p) +Iy, (pq-l) +L, (pea) (5) 
~mx [X,(w-ug+vp) -2Z,(u-vr+wq) ] =M,+M,+M, 


1i+(I,-1,)pq-1,)(P?-q")-I,,(pr +4) +1,(q7-P) 6 
+mx[x,(v+ur-wp) -y (u-vr+wa)]=NytNy+No 
In these equations, the left hand sides represent inertial forces and moments 
(Newton’s Law) and the night hand sides model the external forces. Subscript H reflects 
hydrodynamic contributions, W buoyancy and weight effects, C forces arising from 
control surface (rudders, dive planes, and bow planes) actions, and the rest of the 
symbols are based on standard notation and will be explained at the end of this section. 


The hydrodynamic radiation and viscous forces are expressed as: 


Xy=X,p°+Xq°+X7°+X pr+X +X, wq+X, vpt+X, wr (1) 
+X v7+X w?-Ch ou? 


Yy=V pry s+Y Paty qr+Viv+Y up+yY urty,. vq+Y,,, wp+Y wrt+Y uv+Y vw 


= p [ “Ic p,A(x)(v+xry’ +Cy b(x)(w- ~xq)7] 220) oe ae (8) 
Zy> or p?+Z, pre+Z_y?+Z,, w+Z (9+Z,, vp+Z,. vr+Z,, uw+Z, me 
se aloe h(x)(v+xr)’+Cp b(x)(w- -xq)]-— (w “*Q) ay (9) 

U f* x) 
K,,=K,p+Kr+K gP9 *K,,qr+K,v+K up+Kur+K, vq ill 


+K, wp +K, wr+K uv+K, vw 
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M,=M,q+M_,p?+M,,pr+M,7°+M,w+M_uq+M, vp+M, vr+M,uw+M, v7 


tp fim 24, Sr) ieeeoon Ud (11) 
5° J, Ip, Aa Vv +27)? +Cp be) wag) 1 iar, 


N,=N,ptN F+N pat ,qr+N,v+N up+Nur+N vq+N, w+, wr+N uv+N wr 2) 


1, ft +xr)7+ sop tX!) 
5 Pf IC AGN +27)? +Cp, bea w-xa) 1 aes 


These are given in customary form of series expansion in terms of the 
hydrodynamic coefficients and cross flow integral terms which are integrated over the 


entire length of the body and represent quadratic forces. The cross flow velocity U,, is: 


Uy-[v+zr)? +(w-xq)"]!? (13) 


Hydrodynamic restoring forces and moments are due to the vehicle weight W and 


buoyancy B, and are given by: 


X,,=-(W-B)sin® (14) 
Y,,=(W-B)cos@sind (15) 
Z,y=(W-B)cos@cos} (16) 

Ky=(y.W-y,B)cos8cos@ -(z,, W-z,B)cosOsing (17) 


pad 


My=~-(x,W-x,B)cos6coso -(z-W-z,B)sin8 (18) 


Ny=(x,W-x,B)cosOsind +(yW-y ,B)sin® (19) 


Forces and moments, due to control surface deflections, are reflected as added drag 
in surge, while in sway, heave, pitch, and yaw they are directly proportional to the 


control surface deflections. 


Xo =Uq(X,, 5. +X E 3,0 ) +X,, urd, +X,, uve, +uw(X,, 5, +X,,,5 ») 


q3, s sy ion, 8 2+ Xa 24 mea 6 2) (20) 
Y¥o=¥, ud, (21) 
Z,=uZ,,8,+Z;, 5) (22) 

K,=0 (23) 

M,=u7(M »,5,*M, 5,) (24) 

No=N, u 5 (25) 


Usually, control surface deflections are kept intentionally small, and the linearity 
assumption in Equations (21), (22), (24) , and (25) remains valid. Unlike the surface ship 
case, the roll moment K, is zero for a submersible since the rudder is centered with the 
vehicle center plane. The hydrodynamic coefficients in equations (7) through (12) are 


functions of the frequency of motion, or what amounts to the same thing functions of 


12 





the maneuver at hand. In this work slowly varying reference motions are studied and 
therefore, it is assumed that they remain constant and equal to their zero frequency limit. 
It should be emphasized though, that the constant coefficient assumption would break 
down in studies related to the fast motions under the action of first order wave forces in 
the case of a nearly surfaced submarine. Another important assumption in this study is 
that they are assumed to be constant throughout the range of vehicle angle of attack. 
Ordinary maneuvering models are usually validated for angles of attack + 15 degrees. 
For higher angles of attack, the cross flow drag terms Cp, and Cp, dominate the 
response, and they are functions of the side and angle of attack. Considering them to be 
constant does not alter significantly the behavioral characteristics and qualitative 
bifurcation results that are derived. Similarly , C), and Cp, are functions of speed due 
to the Reynolds number effect on cross flow drag. This is more pronounced in small size 
unmanned untethered vehicles, whereas for submarines the cross flow drag terms remain 
relatively constant over the entire speed range. 

Since the equations of motions are referred to an axis system that is fixed in the 
vehicle, and thus translates and rotates with it, the orientation and the position of the 
moving body axis system relative to fixed inertial reference system must be specified. 
The orientation of the body axis system with respect to the inertial reference system is 
defined by the standard Euler Angles ¥(yaw), O(pitch), and (roll). The rotation 
sequence from the inertial reference system to the body axis system is ¥, 0, ® as shown 
in Figure 4. So to complete the model, we need the expressions for the Euler Angle 


rates of change; 


#3 





(2) Horizontal Flight 
Reference System derived 

from XgYoZqo by rotation a- 
bout Z through yaw angle y. 


(1) Vehicle-Centered 
Gravity-Directed System 
parallel to inertia! 
axis system. 





(3) Roll-Resolved Flight 
Reference System derived 
From X''Y''Z'' by ratation 
about Y'' through pitcaA 
angle 86. 


(4) Vehicle Body Axis Ref- 
erence System derived from 
X'Y'Z' by rotation about X' 


through rol! angle 4. 


Figure 4. Unit Sphere Development of Euler Angles 
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b= p+ q sind tanO + r cosd tan (26) 





6=q cosh-r sind (27) 
v =q sind Pe. cos (28) 
cos8 cos6 


Finally, it is assumed that propulsion 1s inoperative and the propeller is rotating 
freely. Therefore, propulsive forces are not included in Equations (1) through (6). The 
driving mechanism for the vehicle is its excess buoyancy, B-W>Q. The problem is then 
to assess the asymptotic dynamic characteristics of the system during this condition of 
free positive buoyancy ascent. 

The model presented above can be written in its state space form by selecting as 
state variables; 

X,=U ,X%=V ,X,=w , xX, =p (29) 
Xs=Q ,X,=F , =H , x,=0 
where the first six describe the system motion and the last two its geometry. Notice that 
the yaw angle W does not affect the equations of motion and is, therefore, not included 
in (29). Angle Y can be computed from (28), once the time histories of the state 
variables have been obtained. 


In a compact vector form the state equations can be written as, 


x = f(x) | (30) 


where bold face indicates a vector. 
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Major variables and parameters as defined by Smith, Crane, and Summey 


[Ref. 13] are given below: 


1. 


Dynamic Variables 


u,Vv,W - Linear velocity components of vehicle with 
respect to orgin of body axes system relative to fluid. 


p,q,r -Angular velocity components of vehicle with 
respect to body axes system relative to inertial 
reference system. 


X,Y,Z - Hydrodynamic force components along body 
axes. 

K,M,N -Hydrodynamic moment components along body 
axes. 


Mass Distribution Parameters 


m - Mass of the flooded vehicle, including the mass of the 
entrained fluid. 
WwW - Weight of the flooded vehicle , including the weight of 


the entrained fluid. (=gm ;where g is the 
acceleration of gravity). 


V - Displacement volume of the vehicle. 


B - Buoyancy force acting on the vehicle (ogV). This is 
independent of the inertial mass distribution of the 
submersible vehicle, including whether or not it is 
flooded . 


Xo, Yo;2Zq - Coordinates of the CG (center of gravity) in the body 
axis system (Figure 3). These will depend on the mass 
distribution of the vehicle, including the mass of the 
entrained fluid. 


Xp, Yu, Zp - Coordinates of the (center of buoyancy) in the body axis 
system (Figure 3). These are independent of the mass 
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distribution system, but may vary with the addition or 
removal of external appendages. 


| ee re - Moments of inertia about the body system axes, 
including the entrained fluid. 


fF, eel. - Products of inertia about the body system axes, 
including the entrained fluid. 


3. Remaining Parameters 

p - Mass density of fluid medium. 

b(x), h(x) - Width and height of vehicle in its xy and xz 
planes, respectively, at location x measured in the 
body axes system (Figure 3). These quantities are 
required in the integral defining — the cross flow 


forces and moments in the equations of motion. 


Kees Xe - Coordinates of vehicle nose and tail as 
measured in body axis system. (Figure 3). 


6, 64,0, - Sternplane, bowplane and rudder deflection 
angles in radians (Figure 3). 


C. STEADY STATE CONDITIONS 

Steady state conditions are achieved when the submersible reaches constant linear 
and angular velocities which must assume finite values. Therefore, the body fixed linear 
accelerations and body fixed angular accelerations will be zero. Likewise, the vehicle 
will have reached constant angles of roll ¢ and pitch @ , making the derivatives of them 


equal to zero. When these values are put into Equations (1) through (6), and substituting 


6 =6=0 
in Equations (26) through (28) the steady state values of the angular velocities can be 


found as: 
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p =-w sinO (31) 
q = sind cosO (32) 


r = w cos cos 6 (33) 


Combining Equations (31), (32), (33) with the equations (1) through (6) yields a 


system of nine coupled nonlinear algebraic equations in the form 
f(¥) =0 (34) 


where the overbar denotes an equilibrium solution. 

Equation system (34) can be solved for the steady state values of u, v, w, p, g, I, 
@, 9,y. This is a highly nonlinear system of equations, and it may exhibit solution 
branching and/or multiple solutions McKinley [Ref. 14]. Here y is not zero at steady 
State because taking y=O at steady state will restrict the analysis to the vertical plane. 
In this analysis, motion in all six degrees of freedom is studied. 


Equation system (34) can be written in the following form 


fCx,1) = 0 (35) 


Here, denotes any one of the system parameters. By system parameter, we mean the 
effects of control surfaces, propeller rpm the values of buoyancy and weight, etc. 
Therefore, our system of Equations (35) represent a multiparameter problem. This 
multiparameter problem can be solved by keeping all parameters fixed except one 


(denoted by A). This system of equations for various values of can be solved 
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continuation method. The continuation method will be explained in the following 


section. 


D. PRINCIPLES OF CONTINUATION 


The system of nonlinear ’algebraic’ equations 


f(y,4)=0 (36) 


serves as a basis for the discussion. Here y denotes an n-dimensional vector. No 
generality is lost by restricting attention to Equation (36). Both steady state solutions of 
ODE’s and PDE’s are solved by approximating them by such system of nonlinear 
equations. 

Assume that at least one solution of Equation (36) has been calculated. Let us 
denote this first solution by (y' , A,).The continuation problem is to calculate further 


solutions, 


Orn 5 Aas - - - (37) 
until one reaches a target point, say X = A,. (The superscripts are not to be confused 
with exponents.) 
The j-th continuation step starts from (an approximation of) a solution (y’, \,) of 
equation (36) and attempts to calculate the solution (y’*’, \,,,) for the next A, namely 
Nia 


With predictor-corrector methods, the step j>j+1 is split into two steps: 


BO 


(yA) = (97,201) = O71, 4,0) 
Predictor step Corrector step 


see Figure 3 Seydel [Ref. 12]. 





Figure 5. Schematic Showing Predictor Corrector Method (From Seydel, 1988) 


In this, thesis a predictor corrector method is used. In general, the predictor 1s not 
a solution of equation (36). The predictor merely provides an initial guess for corrector 
iterations that home in on a solution of equation (36). The major portion of the work is 
either on the predictor step (resulting in an approximation close to the branch) or on the 
corrector step (if the predictor has produced a guess far from the branch). The distance 
between the two consecutive solutions (y’, ;) and (y’*’, A;,,) is called the step length or 
Step size Seydel [Ref. 12]. 

It is obvious that first of all, we need a solution procedure for the particular class 
of problems under investigation. In this work, a Newton-like method has been used to 


find the solutions of the problem. Therefore, calculation of the first solution is often the 
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most difficult part. Initial guesses may be based on solutions to simplified versions of 
the underlying equations. Sometimes choosing special parameters (often equal to zero) 
leads to equations that are easily solved. In complicated cases, a hierarchy of 
simplifications must be gradually relaxed step by step, the solution of each equation 
serving aS an initial guess to the following more complicated equation. Solving a 
sequence of equations with diminishing degree of simplification until finally the full 
equation solved is called homotopy Seydel [Ref. 12]. 

Equation system (35) can be simplified by assuming the angle of yaw (vy), its 
derivative, and p, q, r to be zero. By doing this, the motion of the vehicle will be 
restricted to the vertical plane. In this case, the resulting equations can be solved 
analytically. These analytical solutions were done in McKinley [Ref. 14]. In this work 
these solutions and certain numerical simulation results have been used as initial 
conditions to the continuation runs. 

The continuation approach does not guarantee that all possible solutions can be 
located in a specific example. In particular, changes to detect isolated branches are 
small. Therefore, in this work there might be other solutions which had not been located 
by the continuation algorithm which is used. In this thesis, a continuation method is 
implemented BIFPACK [Ref. 15]. A careful combination of analytical techniques, 
continuation methods, and numerical simulations are used to ensure that all steady state 


solutions are captured. 


PAL 


Hl. CONTINUATION IN és 


A. GENERAL 

In this chapter, the results of continuation in terms of the sternplane angle will be 
discussed. As was mentioned before, it is very important to have an initial condition 
before starting a continuation study. In his master’s thesis, McKinley solved the system 
of equations in the vertical plane [Ref. 14]. In that work, he found that when the 
longitudinal center of buoyancy (Xgg) was equal to -1 percent of the vehicle length, for 
a certain value of sternplane, deflection the real part of one of the system eigenvalues 
became zero while the others stayed negative. By using this fact, it was assumed that 
there 1s a pitchfork bifurcation in our system for the above specific values of system 
parameters. In this work, the results of his thesis were used as initial conditions in the 
continuation runs to find the vertical plane solutions. Since these results were the 
analytical solutions of the simplified system equations by restricting the motion in vertical 
plane, we were able in this work to validate the continuation runs. In continuation runs, 
the system of equations Equation (35) was solved for changing values of the primary 
parameter A (in this case A=6s). After a few continuation runs, a pitchfork bifurcation 


was found. This pitchfork bifurcation is explained in the following section. 
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B. PITCHFORK BIFURCATION 


Suppose we have the nonlinear system of state equations, 
x = f(x) (38) 
we know that the equilibrium points, x of the system are defined by: 
f (x) =0 (39) 


This is a nonlinear system of algebraic equations, and it may have multiple 
solutions in x, which means that the nonlinear system may have more than one position 
of static equilibrium. If we pick one equilibrium x , we can establish its stability 
properties by linearization. The linearized system becomes 


X=AX : 


where A is the Jacobian Matrix of f(x) evaluated at x, 


*| 


i 
ele 


and the state x has been redefined to designate small deviations from the equilibrium x, 
X7X-X. 


As long as all eigenvalues of A have negative real parts, we know that the linear system 
will be stable. This means that the equilibrium x will be stable for the nonlinear system 
as well. 

The question we ask ourselves next is, what happens if one real eigenvalue of the 


linearized matrix A is zero? The interesting case here is when the rest of the eigenvalues 
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have all negative real parts, otherwise x is unstable and the problem is solved. If the 
case of a zero eigenvalue appears to be too specialized to be of any practical use, 
consider this: Assume that f(x) depends on one physical parameter, and that physical 
parameter is allowed to vary over some range. Then it is clear that A will depend on 
that parameter and as the parameter varies, it is possible that one real eigenvalue of A 
will become zero for a specific value of the parameter. Our problem is then to establish 
the dynamics of the nonlinear system as one real eigenvalue of A crosses zero; 1.e., goes 
from negative to positive. As the solutions evolve with time, things are interesting only 
along the direction of the eigenvector that corresponds to the critical eigenvalue (the one 
that crosses zero). Along the rest of the directions in the state space, everything should 
converge back to equilibrium; remember that we assumed that all remaining eigenvalues 
of A have negative real parts. 

We can see then that it is possible to approximate our original system by a one 
dimensional system, which is much easier to analyze. The dynamics of the two systems 
will be qualitatively similar. The formalization of the above reduction procedure 
constitutes what is known as center manifold reduction, or normal form computation in 
nonlinear analysis. So, let us see what happens for the case of a zero eigenvalue by 


using a (typical) first order system, 
X=Ax-x? , 


for \ > 0 there are two nontrivial equilibria, x= + Vd. The transition of stability is 


illustrated by Figure 6. We can summarize the stability as follows: 
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® For )\ < 0 only the trivial equilibrium exists and is stable. 


@® For) > 0 the trivial equilibrium becomes unstable and a pair of symmetric stable 
equilibria are generated. 







stable 


Stable unstable 


stable 


Figure 6. Supercritical Pitchfork Bifurcation 


This phenomenon, the loss of stability of an equilibrium and the generation of additional 
equilibrium states, is called a pitchfork bifurcation and is very common in nature; Euler 
buckling of a beam is a very typical example. In particular, the above case is referred 
as the supercritical pitchfork, this is a rather benign loss of stability since upon loss of 
stability of the trivial equilibrium the additional nearby equilibrium states are stable. 
Occasionally, the above case is referred to as a soft loss of stability, since for small 
values of \ beyond its critical value, the final steady state of the system does not differ 


much from the nominal (trivial) steady state. 
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As a second example, consider a "similar" system as before, the linear part remains 


the same, and the nonlinear part x* changes sign, 
X=Ax+x? 


equilibria and stability behavior for this equation are shown in Figure 7. There is again 
a loss of stability at the bifurcation point (y, A) = (0, 0), but in contrast to the first 
example of Figure 6 there is no exchange of stability. Instead, the stability is lost locally 
at the bifurcation point. We can summarize the stability as follows : 

® For A > 0 only the trivial equilibrium exists and 1s unstable. 


® For A < 0 the trivial equilibrium becomes stable, and a pair of symmetric 
unstable equilibria are generated. 






=, 


ae 
~~ 
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“~ 
Stable 
— 
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Figure 7. Subcritical Pitchfork Bifurcation 


This case, which is also shown in the figure, is called a subcritical pitchfork. A 


comparison with the first case reveals that this is a much more serious loss of stability 
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case. Upon loss of stability of the trivial equilibrium position, there is no other stable 
equilibrium. 
Now let us perturb the second example, considering the equation 
Ax+x?+B=0 (40) 
An analysis reveals that in the perturbed case the bifurcation is destroyed. 
Differentiating equation (40) with respect to X gives: 


AX 352 AK _ 
x+A a tox Th 0 


Since x=0 is not a solution of equation (40), there is no horizontal tangent 
dx/dA=0 to the solution curve x(A). Checking for a vertical tangent means considering 


the dependence of on x and differentiating with respect to x. This yields 


di 


ae for A=-3x? ’ 


which, substituted equation (40), in tum yields the loci 
pl 2 
Sl eee 
(x,A)=[( 5) ’ Sy ] 
of the vertical tangents. There are three solutions of equation (40) for 
2 
molt Syme 
2 


and the branching diagram looks like Figure 8 (solid lines). The point with a vertical 


tangent is a turning point. For 8 < 0, the curves in Figure 8 are reflected in the A-axis. 


aa 


As 6 — O, the curves approach the pitchfork (dashed curve). For an arbitrary 
perturbation 8, the occurrence of a tuming point is typical. Bifurcation is the rare 


exception because it occurs only for B=0 . 





Figure 8. Perturbed Pitchfork Bifurcation 


In our case, we found a pitchfork bifurcation. We found that the trivial solutions 
of our system come from in-plane (vertical plane) solutions. After the pitchfork 
bifurcation, we had a branching in the steady state solutions. The physical meaning of 
these branches is that they correspond to out of plane motion of the vehicle at the steady 
state. To get these branches, we used the numerical simulation results as an initial guess 
for the continuation run. Figure 9 shows a typical pitchfork bifurcation that we found is 


a supercritical pitchfork bifurcation. In this and all subsequent runs in this chapter, 6s is 
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our main bifurcation parameter (A) and dr is the perturbation parameter (@). In Figure 
9 and all subsequent figures, solid lines correspond to the symmetric system (é6r=0 
degrees) dashed lines to 6r=0.5 degrees, dotted lines to 6dr = 2.5 degrees, and dash-dot 
lines to 6r= 5 degrees. As we can see in this figure, it illustrates all the typical 
characteristics of a pitchfork bifurcation. 


TYPICAL PITCHFORK BIFURCATION FROM CONTINUATION in DS 
0.06 ' 





0.04 
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-20 —15 -—10 =—5 0 +) 10 18 20 


STERN PLANE ANGLE 
Figure 9. Pitchfork Bifurcation From Continuation Runs 


C. CONDITIONS 
1. Defining Additional Terms 


a. Excess Buoyancy 6B 
Excess buoyancy is defined as 6B = B - W where B is the submersible’s 


total buoyancy and W is the submersible’s total weight 
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b. Longitudinal Center of Buoyancy , xc, 
The longitudinal center of buoyancy is defined as xgg = Xg - Xg where 
X¢ is the longitudinal center of gravity with respect to the body fixed axis and x, is the 
longitudinal center of buoyancy with respect to the body fixed axis. 
c. Vertical Center of Buoyancy, Zc, 
The vertical center of buoyancy is defined as Zn = Ze - Z_ Where Z, is 
the vertical center of gravity, with respect to the body fixed axis, and zg is the vertical 
center of buoyancy with respect to the body fixed axis. 


2. Assumed Conditions 


a. Lateral Centers of Gravity, y,, and Buoyancy, y; 
The lateral center of gravity and center of buoyancy are assumed to be 
on the same plane(y, = yz = 0). 
b. Propeller Speed , n (revolution per minute) 
The propeller speed is assumed to be zero (n=O). 
c. Propeller Coefficients , K,,,, and Nprop 


From Smith, Crane, and Summey [Ref. 13], the propeller coefficients are 
ZETO(Kyiop = Norop = 9). 
d. Vertical Center of Buoyancy, Zp 


The vertical center of buoyancy is assumed to be positive. 
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D. STEADY STATE RESULTS 

Figures 10 through 18 show steady state solutions for the eight state variables, 
namely; u, v, w, p, q, r, ®,6. In Figures 10 through 18, solid lines correspond to ér = 
0 degrees, dashed lines correspond to ér =0.5 degree, dotted lines correspond to ér = 
2.5 degrees, and dash-dot lines correspond to dr = 5 degrees. Figures 16 and 17 show 
steady state solutions for roll angle in two different axis scales. In all these cases, 
sternplane action is continuation parameter (i.e. ,A) and dr is the perturbation parameter 
(i.e. ,8). Values of the other parameters are kept fixed:excess buoyancy 6B = 2 % of 
the vehicle weight (W);deflection of bow planes, 5b = 0;location of horizontal and 
vertical centers of buoyancy, x, = Z, = 0; location of vertical center of gravity ,z >, = 


0.1 feet;and location of longitudinal center of buoyancy; xgg= -1 % o the vehicle length. 


SURGE VELOCITY U vs DS 
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Figure 10. Steady State Surge Velocity 
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SWAY VELOCITY V vs OS 
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Figure 11. Steady State Sway Velocity 


HEAVE VELOCITY W vs DS 
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Figure 12. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY p vs DS 
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Figure 13. Steady State Roll Velocity 


PITCH ANGULAR VELOCITY q vs DS 
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Figure 14. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY r vs DS 
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Figure 15. Steady State Yaw Angular Velocity 


ROLL ANGLE PHI vs DS 


PHI 





—720 ~ 5 =e ee) 0 a 10 15 20 
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Figure 16. Steady State Roll Angle 
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ROLL ANGLE PHI vs DS 
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STERN PLANE ANGLE 
Figure 17. Steady State Roll Angle 


PITCH ANGLE THETA vs DS 
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Figure 18. Steady State Pitch Angle 
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It can be seen that all horizontal plane variables, v, p, r and ®, exhibit the typical 
characteristics of pitchfork bifurcation. For ér= 0, there exists a dive plane angle ds 
between -5 and -10 degrees, such that the in plane solution becomes unstable and 
symmetric stable out of plane solutions appear. Nonzero values of dr result, naturally, 
in out of plane solutions for the entire range of ds. Multiple solutions appear also at a 
certain value of 5s, called turning point, which is a function of 6dr. 

The vertical plane variables u, w, q, and 6 are even functions of horizontal plane 
motions and, as a result, they exhibit "solution separation" of the in plane motions for 
values of 5s lower than the bifurcation point. Of special interest is an examination of the 
steady state values of the pitch angle @ and roll angle ®, since these are directly related 
to the vehicle orientation. A clear examination of the possible steady state solution, so 
is the pair (z-6, a-%). When dr= Q, it can be seen from Figure 18 that the stable in- 
plane solution exhibits a steady state pitch angle greater than 90 degrees for values of ds 
less than about -5 degrees. This is, of course, true up to the pitchfork bifurcation point. 
This phenomenon of inverted pendulum stabilization was first discussed by McKinley 
[Ref.14]. It appears, however, from the results of Figure 18 that the inverted pendulum 
stabilization is highly degenerate and is destroyed as soon as a non-zero rudder angle ér 
is applied. The two solutions, 6 and z-6, do not cross each other as they did for ér=0. 
Instead, they veer of each other, the lower one corresponding to =O and the upper one 
to =z. Asa result, the vehicle configuration is always keel down, the only difference 


between the two solution sets at Figures 17 and 18 being 180 degrees heading. 
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Another important observation can be made by considering the results of Figures 
16 and 17, in terms of the roll angle ®. It appears that there exists a range of sternplane 
angle where the exact value of ® is very sensitive to the exact value of é6r for small 
values of ér. The rudder angle is related to the angle of attack of the fluid flow, which 
is highly dependent on such factors as vortex shedding, flow separation, and ambient 
flow turbulence Sarpkaya [Ref. 5]. Since these factors are to a great extent uncertain and 
difficult to model, situations like the one described above should be avoided in practice. 
This bifurcation and continuation methods used in this thesis provide a systematic and 
effective way of identifying precisely parameter regions where such sensitivity of the 
response to parameter values exists. These regions should be avoided when executing a 


proper recovery procedure. 
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IV. CONTINUATION IN or 


A. GENERAL 

In this chapter, the results of continuation in terms of the rudder angle will be 
discussed. In this continuation run, the initial conditions were taken from the previous 
results of continuation in terms of the rudder angle. The configuration of the vehicle was 
just like the previous case. In this case the system of equations (35) was solved for 
changing values of the primary parameter A (here A= dr). At the end of the runs, a 
hysteresis behavior was observed for some of the state variables. This hysteresis 
behavior, which is complementary to pitchfork bifurcation, is explained in the following 


section. 


B. HYSTERESIS PHENOMENON 
When one-parameter families of equation are studied, there is no need to talk about 
a hysteresis phenomenon. But similarly to the pitchfork bifurcation, this situation might 


Change when a two-parameter family of equations is studied, 
f(y,a,a@) 


Now let us consider the following equation: 
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x3-h + ax = 0 (41) 


Let us call this algebraic equation as a function of x F(x) = Q. It can be observed that 
equation (41) is similar to the generic pitchfork equation (40). The only difference is that 
the roles of } and a have been interchanged. The constant term is now our primary 
bifurcation parameter, while the coefficient of x is the perturbation parameter. We know 


that at the turning point F, and dF/dx, should be zero. So, 


ead — 0) ==> xe=- 


wile 


If we substitute this value into Equation (41), we will get the following equation. 


27A7+4a°=0 (42) 


This equation gives the relationship between ) and a at the vertical tangent. Now we 


Figure 19. Cusp Curve 


consider equation (41) as a two parameter model with A and a having equal nights. 
Equation (42) establishes a curve in the (a , A) plane, which takes the form of a cusp 
(Figure 19), (Cusp curves will be explained in a more general context in Chapter VII.). 

For all combinations of ( a, A) values inside the cusp (hatched region), Equation 
(41) has three solutions; for values outside the cusp, there is only one solution. This 
becomes clear when branching diagrams are drawn. As a generic behavior, if we draw 
branching diagrams for different values of @ when A is the primary continuation 


parameter, we get the Figures 20, 21, 22 according to the sign of a. 


®a>d0 


®a<0 





Figure 20. Branching Diagram When a@ = 0 
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Figure 21. Branching Diagram When a < 0 


Figure 22. Branching Diagram When a > 0 


In this last case, we see a hysteresis effect. When hysteresis occurs in a branch, 
there are two associated turning points with hysteresis behavior. These points occur 
precisely at the two real solutions of \ from Equation 42. Here there is a jump when the 


solution is on an appropriate part of the branch. The two outer solution in x are stable 
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while the inner most solution between the two turning points is unstable. The two dotted 
straight lines of Figure 19 represent two distinct paths through the cusp and explained 


in Chapter VHT. 


C. STEADY STATE RESULTS 

Figures 23 through 30 show steady state solutions for the eight state variables, 
namely; u, v ,w, p, q, r, ®, 8. In Figures 23 through 30 solid lines correspond to 6s= 
0 degree, dashed lines to 6s = -5 degrees, dotted lines to 6s = -10 degrees and dash-dot 
lines to 6s = -20 degrees. For these continuation runs, the rudder angle (6r) is the 
primary continuation parameter (i.e,A) and stern plane angle (6s) 1s perturbation 
parameter (i.e,a) and values of all other parameters are kept fixed:excess buoyancy, 
6B = 2 % of the vehicle weight (W); deflection of bow planes, 6b = 0 ; location of 
horizontal and vertical centers of buoyancy ,x,= Z, = O;location of longitudinal center 


of buoyancy, Xcp =~! % of the vehicle length (L). 
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SURGE VELOCITY U vs DR 
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Figure 23. Steady State Surge Velocity 


SWAY VELOCITY V vs DR 
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Figure 24. Steady State Sway Velocity 
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HEAVE VELOCITY W vs OR 
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Figure 25. Steady State Heave Velocity 


ROLL ANGULAR VELOCITY P vs OR 
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Figure 26. Steady State Roll Angular Velocity 
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Figure 27. 
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Figure 28. 


PITCH ANGULAR VELOCITY U vs DR 
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Steady State Pitch Angular Velocity 


YAW ANGULAR VELOCITY R vs DR 
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Steady State Yaw Angular Velocity 
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ROLL ANGLE PHI vs DR 





Figure 29. 


THETA 


Figure 30. 
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Steady State Roll Angle 


PITCH ANGLE THETA vs OR 
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Steady State Pitch Angle 
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Figures 23 through 30 show that the horizontal plane state vanables v, p, r, ® 
behave like Figure 22 and exhibit a hysteresis behavior for certain values of 6s. If we 
zoom in Figure 28, we obtain Figure 31 and can see this hysteresis effect very clearly. 
When és = -15 degrees, there is a typical hysteresis phenomenon. Figures 24, 26, and 


28 also show the same kind of behavior. 


TYPICAL HYSTERESIS BEHAVIOR FROM CONTINUATION RUNS 


-0.04 


-0.06 
= 





RUDDER ANGLE 
Figure 31. Typical Hysteresis Effect from Continuation 


On the other hand, the vertical plane state variables u, w, q, @ behave like an even 
function in x (for example | x | ). This behavior can be seen in Figures 23, 25, 27, 30. 
where it is clear that they exhibit generic behavior of an even function for different 
values of the perturbation parameter (i.e, 6s). If we zoom in Figure 23, we can see this 


behavior, very clearly. 


4‘ 


TYPICAL BEHAVIOR OF AN EVEN FUCTION FROM CONTINUATION 





RUODER ANGLE 
Figure 32. An Even Function from Continuation 
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V. CONTINUATION IN xc, 


A. GENERAL 

In this chapter, the results of continuation in terms of the longitudinal center of 
buoyancy will be discussed. In continuation runs, the system of equations (35) was 
solved for changing values of the primary parameter \ where in this case A = Xgp,). In 
these runs, first we fixed the value of rudder angle and changed the values of stern plane 
angle, then we kept the sternplane angle fixed and changed the rudder angle, to see how 
the steady state values were changing. When we keep the rudder angle at zero and 
change the values of sternplane angle, we observe that there is a pitchfork bifurcation for 
the states v, p, r. Then, we fixed ds at -20 degrees and started to change rudder angle 
to see how the pitchfork bifurcation diagrams were changing with changing rudder angle. 
As a last continuation run for this configuration, the value of sternplane angle was fixed 
at zero and we changed the value of rudder angle. We didn’t see any solution branching 


for this case. 


B. STEADY STATE RESULTS 


1. 6dr=0 Degrees, és is Perturbation Parameter 
Figures 33 through 40 show steady state solutions for eight state variables, 
namely; u, v, W, p, q, r, ®, 6. In Figures 33 through 40, solid lines correspond to 6s=0 


degrees, dotted lines correspond to 6s=-5 degrees, and dash-dot lines correspond to és=- 
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20 degrees. In graphs of the state, variables v, p, r, 6 dash-dot lines may seem like solid 
lines; this is because the continuation between the parameter values xg, = -13 degrees 
and Xo, = .7 degrees traced the same curve twice in order to make sure that all solutions 
for all the state variables were captured. Therefore, out of plane solutions in this 
parameter range belong to sternplane angle 6s = -20 degrees. 

In all these cases, the longitudinal center of buoyancy is the continuation parameter 
and 6s is the perturbation parameter, while dr is kept at zero degrees. Values of the other 
parameters are kept fixed: excess buoyancy, 6B =2 % of the vehicle weight (W); 
deflection of the bow planes, 6b = Q: location of horizontal and vertical centers of 


gravity, Xx, = Z, =0: location of vertical center of gravity, z4,= .1 feet. 


SURGE VELOCITY U vs XGB 


eeeorereeoeeeoe dP eC eoeeeseeeeeemPee ese eeeeeaate ewe eo + oe eo 





e Med rem 
-” cam ° 
° - a a = 
s = 
“om «a? 
ie 


ae 
ae -1.5 —1 —0.5 0 eR) 1 beg 2 


LONGITUDINAL CENTER OF BUOYANCY %L 
Figure 33. Steady State Surge Velocity 
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SWAY VELOCITY V vs XGB 
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Figure 35. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs XGB 
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Figure 36. Steady State Roll Angular Velocity 
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Figure 37. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY R vs XGB 
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Figure 38. Steady State Yaw Angular Velocity 


ROLL ANGULAR VELOCITY PHI vs XGB 
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Figure 39. Steady State Roll Angle 
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PITCH ANGULAR VELOCITY THETA vs XGB 
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Figure 40. Steady State Pitch Angle 


2.  6s=-20 Degrees, dr is the Perturbation Parameter 

Figures 41 through 48 show steady state solutions for the state variables u, 
v, W, p, q, r, ®, 8. In Figures 41 through 48, solid lines correspond to ér=O degrees, 
dotted lines correspond to 6r=2.5 degrees, and dash-dot lines correspond to ér=5 
degrees. 

In all of these cases, longitudinal center of buoyancy is the primary 
continuation parameter and ér is the perturbation parameter, while és is kept at -20 
degrees. Values of the other parameters are kept fixed: Excess buoyancy, 6B = 2 % of 
the vehicle weight (W): deflection of bow planes, 6b = 0; location of horizontal and 


vertical centers of buoyancy, x, = Z, =0; location of the vertical center of gravity, z4= 


.1 feet. 
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SURGE VELOCITY U vs XGB 
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Figure 41. Steady State Surge Velocity 


SURGE VELOCITY V vs XGB 
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Figure 42. Steady State Sway Velocity 
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HEAVE VELOCITY W vs XGB 
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Figure 43. Steady State Heave Velocity 
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Figure 44. Steady State Roll Angular Velocity 
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PITCH ANGULAR VELOCITY @ vs XGB 
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Figure 45. Steady State Yaw Angular Velocity 


YAW ANGULAR VELOCITY R vs XGB 
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Figure 46. Steady State Yaw Angular Velocity 


37 





YAW ANGLE PHI vs XGB 
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Figure 47. Steady State Roll Angle 


PITCH ANGLE THETA vs XGB8 
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Figure 48. Steady State Pitch Angle 
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3. 6s = 0 Degrees, dr is the Perturbation Parameter 

Figures 49 through 56 show steady state solutions for the eight state 
variables, namely u, v, w, p, q, r, ®, 8. In Figures 49 through 56, solid lines 
correspond toédr = 0 degrees, dashed lines correspond to ér = 0.5 degrees, dotted lines 
correspond to 6dr = 2.5 degrees and dash-dot lines correspond to ér = 5 degrees. 

In all these cases, the longitudinal center of buoyancy is the primary 
continuation parameter and 6dr is the perturbation parameter. Values of the other 
parameters are kept fixed: Sternplane action, 6s = 0 degrees:excess buoyancy, 6B = 2 
% of the vehicle weight (W); deflection of bow planes, 6b = 0; location of horizontal 
and vertical centers of buoyancy, x, = Zg, = 0; location of vertical center buoyancy, Zc, 


= 0.1 feet. 
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Figure 49. Steady State Surge Velocity 
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SWAY VELOCITY V vs XGB 
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Figure 50. Steady State Sway Velocity 


HEAVE VELOCITY W vs XGB 
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Figure 51. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs XGB 
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Figure 52. Steady State Roll Angular Velocity 
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Figure 53. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY R vs XGB 
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Figure 54. Steady State Yaw Angular Velocity 
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Figure 55. Steady State Roll Angle 
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PITCH ANGLE THETA vs XG8 
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Figure 56. Steady State Pitch Angle 
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VI. CONTINUATION IN 6B 


A. GENERAL 

In this chapter, the results of continuation in terms of the excess buoyancy will be 
discussed. In continuation runs, the system of equations (35) was solved for changing 
values of the primary parameter ) (in this case A= 6B). In these runs, first we fixed the 
value of rudder angle and changed the values of sternplane angle, then we kept the 
sternplane angle fixed at zero and changed the rudder angle to see how the steady state 
values were changing. When we kept the rudder angle at zero and changed the values of 
sternplane angle, we saw that there was a pitchfork bifurcation for the state variables v, 
p, r, and ®. After we observed this pitchfork bifurcation for 6s= -20 degrees and ér=0 
degrees, we started to change the rudder angle dr to see how the pitchfork bifurcation 
diagrams were changing. As an other continuation run we fixed ér=0 degrees and 
changed 6s towards positive values; however, we didn’t see any branching in this case. 
As a last continuation run for this configuration, the value of sternplane angle was fixed 
at zero and the value of rudder angle was changed. We didn’t see any branching behavior 


for this case. 


B. STEADY STATE DIAGRAMS 


1. dr=0 Degrees, ds is the Perturbation Parameter 
Figures 57 through 64 show steady state solutions for eight state variables, 
namely; u, v, w, p, q, r, ®, @. In Figures 57 through 64, solid lines correspond to 6s=0 
degrees, dashed lines correspond to 6s=-5 degrees, and dash-dot lines correspond to 
6s=-20 degrees. In all of these cases, excess buoyancy is the continuation parameter, 
while é6r is kept at zero degrees. Values of the other parameters are kept fixed: 
longitudinal center of buoyancy x,,=-1 % of the vehicles length (L); deflection of the 
bow planes, 6b=0; location of horizontal and vertical centers of buoyancy, x,=z,=0; 
location of vertical center of buoyancy, z,,=.1 feet. 
As we can see from Figures 57 through 64, there is a pitchfork bifurcation for state 
variables u, v, p, r, ®. We have out of plane solutions for sternplane angle 6s=-20 


degrees, this is also the configuration in which we observe bifurcations. 
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SURGE VELOCITY U vs DELB 
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Figure 57. Steady State Surge Velocity 
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SWAY VELOCITY V vs DELB 
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Figure 58. Steady State Sway Velocity 
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Figure 59. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs DELB 
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Figure 60. Steady State Roll Angular Velocity 


PITCH ANGULAR VELOCITY Q vs DELB 
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Figure 61. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY R vs DELB 
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Figure 62. Steady State Yaw Angular Velocity 
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Figure 63. Steady State Roll Angle 
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PITCH ANGLE THETA vs DELB 
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Figure 64. Steady State Pitch Angle 


After analyzing the above diagrams, we decided to look at a case with a 
positive stemplane angle. We observed that there was no solution branching for this 
case. We got the inplane solutions and there were no out of plane solutions. In Figure 
65, we can see the situation for the state variable u. In this figure, solid lines correspond 
to 6s=0 degrees, dotted lines correspond to 6s=5 degrees, and dash-dot lines correspond 
to ds= 20 degrees. The steady state solutions are not very sensitive to positive sternplane 


action. 
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SURGE VELOCITY U vs DELB 
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Figure 65. Steady State Surge Velocity When 6s is Positive 


2. 6s = -20 Degrees, dr is the Perturbation Parameter 

Figures 66 through 73 show steady state solutions for the state variables u, 
v, W, p, q, r, ®, @. In these figures, dashed lines correspond to é6r= 0.5 degrees, and 
dash-dot lines correspond to ér= 5 degrees. In all of these the cases, excess buoyancy 
is the primary continuation parameter and dr is the perturbation parameter, while ds is 
kept -20 degrees. Values of the other parameters are kept fixed. Longitudinal center of 
buoyancy, -1% of the vehicle length (L); deflection of bow planes, 5b= 0; location of 
horizontal and vertical centers of buoyancy, x,= Z,= OQ; location of the vertical center 
of buoyancy, Zgg= 0.1 feet. By doing this, we saw how the bifurcation points were 


changing as we changed the perturbation parameter 6dr. 
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Figure 66. 
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Steady State Surge Velocity 


SWAY VELOCITY V vs DELB 
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Figure 68. Steady State Heave Velocity 
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Figure 69. Steady State Roll Angular Velocity 
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PITCH ANGULAR VELOCITY Q vs DELB 
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Figure 70. Steady State Pitch Angular Velocity 
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Figure 71. Steady State Yaw Angular Velocity 
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ROLL ANGLE PHI vs DELB 
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Figure 72. Steady State Roll Angle 
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Figure 73. Steady State Pitch Angle 
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3.  6s=0 Degrees, dor is Perturbation Parameter 

Figures 74 through 81 show steady state solutions for the eight state 
variables, namely u, v, w, p, q, r, ®, 6. In Figures 74 through 81 solid lines correspond 
to 6dr= O degrees, dashed lines correspond to 6r= 0.5 degrees, and dash-dot lines 
correspond to ér=5 degrees. In all these cases, the excess buoyancy 6B is the primary 
continuation parameter and or is the perturbation parameter. Values of the other 
parameters are kept fixed: Sternplane action, 6s=0 degrees; longitudinal center of 
buoyancy xgg=-1% of the vehicle length (L); deflection of bowplanes, 6b=0; location 
of horizontal and vertical centers of buoyancy, x,= Z,= 0; location of vertical center of 
buoyancy, Zgg= 0.1 feet. As we can see, there was no solution branching for this case. 


Of course, we got out of plane solutions for nonzero rudder angle values. 
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Figure 74. Steady State Surge Velocity 
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Figure 75. 
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Figure 77. Steady State Pitch Angular Velocity 
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Figure 78. Steady State Roll Angle 
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YAW ANGULAR VELOCITY R vs DELB 
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Figure 79. Steady State Yaw Angular Velocity 
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Figure 80. Steady State Roll Angle 
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PITCH ANGLE THETA vs DELB 
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Figure 81. Steady State Pitch Angle 
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VI. CONTINUATION IN 2, 


A. GENERAL 

In this chapter, the results of continuation in terms of the vertical center of 
buoyancy will be discussed. In continuation runs, the system of equations (35) was solved 
for changing values of the primary parameter A (in this case A= Zg,). In these runs, first 
we fixed the value of rudder angle and changed the values of sternplane angle, then we 
kept the sternplane angle fixed at zero and changed the rudder angle to see how the 
steady state values were changing. When we kept the rudder angle at zero and changed 
the values of sternplane angle, we saw that there was a pitchfork bifurcation for the state 
variables u, v, w, p, r, ®, 6. After we saw that there was a pitchfork bifurcation when 
ds= -20 degrees and ér= 0 degrees, we started to change rudder angle dr to see how the 
pitchfork bifurcation diagrams were changing. As a last continuation run for this 
configuration, the value of stern plane angle was fixed at zero and the value of rudder 


angle was changed. We didn’t observe any branching behavior for this case. 


B. STEADY STATE RESULTS 


1. 6d6r= 0 Degrees, és is the Perturbation Parameter 
Figures 82 through 89 show steady state solutions for eight state variables, 
namely; u, v, w, p, q, r, ®, @. In Figures 82 through 89, solid lines correspond to 6s= 


O degrees, dashed lines correspond to 6s=-5 degrees and dash-dot lines correspond to 
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6s=-20 degrees. In all of these cases, the vertical center of buoyancy is the primary 
continuation parameter, and és is the perturbation parameter, while dr is kept at zero 
degrees. Values of the other parameters are kept fixed: Excess buoyancy 6B=2% of the 
vehicle weight (W); longitudinal center of buoyancy x,,=-1% of vehicle length; 
deflection of the bow planes 5b=0; location of the vertical and horizontal centers of 
buoyancy, Zg=Xg= 0. As we see from Figures 82 through 89, there is a pitchfork 
bifurcation for state variables u, v, w, p, r, ®, 9. We have out of plane solutions for 


this vehicle configuration. 
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Figure 82. Steady State Surge Velocity 
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SWAY VELOCITY V vs ZGB 
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Figure 83. Steady State Sway Velocity 
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Figure 84. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs ZGB 
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Figure 85. Steady State Roll Angular Velocity 
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Figure 86. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY R vs ZGB 
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Figure 89. Steady State Pitch Angle 


2.  6s=-20 Degrees, dr is the Perturbation Parameter 

Figures 90 through 97 show steady state solutions for the state variables u, 
Vv, W, p, q, fT, &, 6. In Figures 90 through 97, solid lines correspond to é6r=0 degrees, 
dashed lines when 6r=0.5 degrees, dash-dot lines correspond to dr=5 degrees. In all 
of these cases, vertical center of buoyancy Zg, is primary continuation parameter and dér 
is the perturbation parameter, while 5s is kept at -20 degrees. Values of the other 
parameters are kept fixed; excess buoyancy 6,=2% vehicle weight (W), longitudinal 
center of buoyancy xg,=-1% of the vehicle length (L); deflection of bowplanes 6b=0; 
location of horizontal and vertical centers of buoyancy, x,=z,;=0. By doing this, we 


saw how the bifurcation points were changing as we change perturbation parameter or. 
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A bifurcation phenomenon took place for 6r=0 degrees and disappeared after we applied 


a small nonzero rudder angle value. 


SURGE VELOCITY U vs ZGB 





0 0.05 0.1 Oo Or Orzo OES 


VERTICAL CENTER OF BUOYANCY (feet) 
Figure 90. Steady State Surge Velocity 
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Figure 91. Steady State Sway Velocity 


87 


HEAVE VELOCITY W vs ZGB 
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Figure 92. Steady State Heave Velocity 


ROLL ANGULAR VELOCITY P vs ZGB 


0.06 


0.02 


-0.04 





-0.06 
fe) 0.05 Ort 0.15 0.2 0.25 0.35 


VERTICAL CENTER OF BUOYANCY (feet) 
Figure 93. Steady State Roll Angular Velocity 
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PITCH ANGULAR VELOCITY Q vs ZGB 
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Figure 94. Steady State Pitch Angular Velocity 
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Figure 95. Steady State Yaw Angular Velocity 
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ROLL ANGLE PH! vs ZGB 
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Figure 96. Steady State Roll Angle 
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Figure 97. Steady State Pitch Angle 
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3. 6s = 0 Degrees, dr is the Perturbation Parameter 

Figures 98 through 105 show steady state solutions for the eight state 
variables, namely u, v, w, p, q, r, ®, 9. In Figures 98 through 105, solid lines 
correspond to 6r=0 degrees, dashed lines correspond to 6dr=.5 degrees, dash-dot lines 
correspond to 6dr=5 degrees. In all these cases, the vertical center of buoyancy is the 
primary continuation parameter. Values of the other parameters are kept fixed; sternplane 
action 6s=0 degrees; longitudinal center of buoyancy xg, = -1% of vehicle length (L); 
deflection of mow planes, 5b=0 degrees; excess buoyancy B= 2% of vehicle weight 


(W); location of horizontal and vertical centers of buoyancy x,=z,=0. As we see, there 
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Figure 98. Steady State Surge Velocity 
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is no solution branching for this last case. We have out of plane solutions for nonzero 


rudder angles. 
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Figure 99. Steady State Sway Velocity 
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Figure 100. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs ZGB 
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Figure 101. Steady State Roll Angular Velocity 
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Figure 102. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY R vs ZGB 
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Figure 103. Steady State Yaw Angular Velocity 
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Figure 104. Steady State Roll Angle 
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PITCH ANGLE THETA vs ZGB 
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Figure 105. Steady State Pitch Angle 
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Vil. CUSP CURVES 


A. GENERAL 

In this chapter, the cusp phenomenon will be discussed. After seeing the branching 
characteristics of the steady state solutions of our system, we decided to look at how the 
bifurcation points and the turning points were changing with changing values of rudder 
angle. The main difference between a bifurcation point and a tuming point is that, at a 
bifurcation point, exactly two branches intersect with two distinct tangents; whereas, at 
a turning point, the branch comes from one side and turns back. We can see this 
phenomenon from Figure 106. This is the graph that was generated for the surge velocity 
when the prmary continuation parameter was X,g, and the sternplane angle 6s=-20 
degrees. In Figure 106, points A and B are bifurcation points, and points C and D are 
turning points. Clearly, locally there are no solutions on one side of a turning point and 
two solutions on the other side, but there are solutions on each side of a bifurcation 
point. 


Let us look at Equation 43 again, this is the model equation that represents a cusp, 
x?-Bx+A=0 (43) 


Let us call this algebraic equation as a function of x, F(x)=0. For the critical curve, F 


and dF/dx should be zero. 
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TURNING AND BIFURCATION POINTS 


SURGE VELOCITY 





Figure 106. Turning and Bifurcation Points 


CUSP CURVE 


Figure 107. Cusp Curve 
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If we substitute this value into equation (43), we get the following relation: 
4B?=27A? 


If we plot this relation in the AB plane, we get the cusp curve shown in Figure 107. This 
represents a separation. Equation 43 has three solutions for parameter values inside the 
cusp and one solution outside this parameter values. In Figure 107, if we fix the value 
of A and vary the value of B (path (a) in Figure 107) we get a pitchfork bifurcation 


diagram, as shown in Figure 108. 


CONSTANT A, CHANGING 8 
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Figure 108. Pitchfork Bifurcation 
In Figure 107, if we fix the value of the value of B and vary the value of A (path (b) in 


Figure 107) we get a hysteresis behavior, as shown in Figure 110. It is clear, that when 
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a line crosses the cusp curve in the (A, B) plane a branching phenomenon takes plane. 


As mentioned before, Equation 43 has three solutions for parameter values inside the 


cusp and one solution outside the cusp this can be seen easily in Figure 109. Now let us 


take a look at Figure 111. As can be seen, there is a simple way of recovering a 


bifurcation diagram from the path through the given cusp by lifting the path to the cusp 


surface. 
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Figure 110. 


Figure 111. 
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Geometry of Cusp Phenomenon 
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If we have the cusp curves for a multiparameter equations system, we can always 
have an idea about the solutions of this system for different perturbed cases. This means 
that by looking at the path of the parameters through the cusp curve, we can see how 
many solutions we have for the corresponding configuration of the system. Similarly, we 
can think about another problem which has a turning point. If we draw the position of 
turning point in the plane with changing parameter values, we will get two solutions 
inside the curve and no solutions outside these parameter values. Therefore, these graphs 
can tell us the nature of the solution set of our system. In our study, we tried to get the 
cusp curves and the position of the tuming points and the results are presented in the 


following section. 


B. RESULTS 

First of all, let us look at the case when we have the sternplane action as our 
primary continuation parameter. When ds is between -5 and -10 degrees, we had a 
bifurcation point and in that case the rudder angle was zero, see Figure 13. To get the 
complete cusp curve, we changed the value of dr and found the locations of the 


bifurcation point. These were plotted as ds versus dr, see Figure 112. 
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CUSP CURVE WHEN CONTINUATION PARAMETER is OS 
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Figure 112. Cusp Curve When Continuation Parameter is 
Sternplane angle 


In Chapter V, we saw that when we have xg, aS primary continuation parameter, 
and sternplane action é6s=-20 degrees, we have both bifurcation and tuming points, 
Figure 106. To see the effect of changing the rudder angle on the bifurcation and tuming 
points, we changed the value of rudder angle and found the location of bifurcation and 
turning points in the (xg, dr) plane. The results are shown on the following graphs. The 
bifurcation points generated a cusp curve and the tuming points generated a parabola in 


the (Xg,,6r) plane. 
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Figure 113. Cusp Curve When Continuation Parameter iS Xg, 


TURNING POINT 
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Figure 114. Turning Point When Continuation Parameter is xg, 
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CUSP CURVE 





=O fe) Cal 0.2 0.5 0.4 0:5 0.6 Oa7 


XGB 
Figure 115. Second Cusp Curve When Continuation Parameter is 
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Figure 116. Second Turning Point When Continuation Parameter 
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In Chapter VI, we used the excess buoyancy as our primary continuation 
parameter. In that configuration, we had a tuming point and a bifurcation point for 
sternplane angle 5s=-20 degrees, (Figure 57 in Chapter VI). To see the effect of 
different rudder angles on these points, we changed the value of the rudder angle and 
found the location of these points in the (Xg,, dr) plane. These results are shown in 


Figure 117 and Figure 118. Again, we got a cusp curve and a parabola. 
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Figure 117. Cusp Curve When Continuation Parameter Is Excess 
Buoyancy 


105 


TURNING POINT 
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Figure 118. Turning Point When Continuation Parameter is 


Excess Buoyancy 


In Chapter VII, we used the vertical center of buoyancy as our primary 
continuation parameter. In that continuation, run we found that there were two 
bifurcation points at sternplane angle 6s=-20 degrees, see Figure 82 Chapter VIL. In that 
figure, the bifurcation point when Z,, is around zero is not clear, but it becomes clear 
when we changed the value of rudder angle. To construct the cusp curves, we changed 
the value of rudder angle and found the coordinates of bifurcation points in the (Zp, 6r) 


plane. These graphs are shown in Figures 119 and 120. 
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CUSP CURVE 
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As a final study in this chapter, we changed the value of y, in the first continuation 
run and tried to get the biased cusp curves for this perturbed case. As a result, we can 
still see a symmetric bifurcation for this perturbed case for a specific parameter 
combination. It is clear, however, that this situation is very degenerate and can easily be 
destroyed for different parameter values. In Figure 121, solid lines correspond to yg=0 
feet, dashed lines correspond to y,=0.001 feet, dotted lines correspond to y,=0.003 


feet, and dash-dot lines correspond to yg=0.005 feet. 
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Figure 121. Cusp Curve When Continuation Parameter is 6s for 
Different Values of y, 
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IX. DYNAMIC STABILITY ANALYSIS 


A. GENERAL 

So far, we have performed continuation runs and have obtained steady state 
solutions of the system. For stability analysis, we chose a representative case and 
performed the stability analysis for that case. We selected the case when we had xg, as 
Our primary continuation parameter and rudder angle 6r=0 degrees, and sternplane angle 
6s=-20 degrees. To predict dynamic stability, we used the six degree of freedom 
equations of motion along with the Euler angle rate equations for the derivatives of the 
angles of pitch and roll (@ and @). These equations of motion were then linearized around 
the steady state nominal points computed in the previous chapters. Eigenvalue analysis 


was then used to compute the stability characteristics of each solution. 


B. DYNAMIC STABILITY RESULTS 

The stability analysis results are shown in Figure 122 and Figure 123. In these 
figures, solid lines correspond to stable solutions while dashed lines correspond to 
unstable solutions. In Figure 122 up to point A, we have a stable solution, which is 
predominantly forward motion, because we have a high surge velocity and the vehicle 
is moving in the vertical plane. Between point A and B an out of plane solution exists 


and the vehicle exhibits motion in the horizontal plane, also. From Figure 123, it is clear 
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that at point A the sway velocity becomes nonzero. This out of plane solution is stable 


up to point B, where 


SURGE VELOCITY 





Figure 122. Steady State Surge Velocity 


SWAY VELOCITY 





Figure 123. Steady State Sway Velocity 
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the out of plane motion becomes unstable up to point C, where it disappears. In Figure 
122, the motion of the vehicle from point D to E is stable and has a low surge velocity, 
which means that the motion of the vehicle is predominantly in the upward direction. 
The vehicle is moving in the vertical plane and there is no horizontal plane motion. 
To assess stability of each solution, we used eigenvalue analysis. The eigenvalues 
that we computed for the following cases are summarized in Table 1. 
@ Case | 
Sternplane angle 6s=-20 degrees, rudder angle 6r=O degrees, xgg=-2 % L, in 
plane solution. 
@ Case 2 
Sternplane angle 6s=-20 degrees, rudder angle 6r=0 degrees, x,,=-1.25 % L, out 
of plane solution. 
® Case 3 
Sternplane angle 6s=-20 degrees, rudder angle r=O degrees, x,,=0.57% L point 
C Figure 122. 
@ Case 4 
Sternplane angle 6s=-20 degrees, rudder angle dr=O degrees, x,,=0 % L. This 
is on the stable branch between point D and E in Figure 122. 
Table 1 represents the eigenvalues of the system. We can see for case three that 
there is a pair of complex conjugate eigenvalues with positive real parts. This means that 
we have a hopf bifurcation prior to that point, and we expect periodic solutions for that 


part of the solution branch. Since we have a stable solution up to point E in Figure 122, 
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we don’t expect to see this periodic out of plane behavior up to point E. The methods 
that we used in this work give only stationary solutions, therefore, to see this periodic 
behavior we used numerical simulation techniques. 


TABLE 1 


-0.198-0.115 +0.015-0.004 














C. NUMERICAL INTEGRATION RESULTS 


The linearized dynamic response results were verified by simulations using 
numerical integration of the full six degrees of freedom equations of motion for the 
swimmer delivery vehicle (SDV). Figures 124 and 125 show a plot of surge velocity 
(U), and sway velocity (V) versus time respectively for the center of gravity aft of the 


center of buoyancy case (xg,=-2) with a diveplane angle (és) of -20 degrees, and rudder 
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Figure 124. Time History of Surge Velocity for Case 1. 
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Figure 125. Time History of Sway Velocity for Case 1. 
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angle (6r) of 0 degrees. The steady state results of the numerical integration method 
match the linearized dynamic results exactly. 

Figures 126 and 127 show a plot of surge velocity (U), and sway velocity (V) 
versus time respectively for the center of gravity aft of the center of buoyancy case 
(Xg6g= -1.25) with a diveplane angle (6s) of -20 degrees and rudder angle (dr) of 0 
degrees. The steady state results of the numerical integration method match the linearized 
dynamic results exactly. Vehicle has an out of plane motion 1n this case. As can be seen 


in Figure 126. the sway velocity V converged to a nonzero Steady state value. 
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Figure 126. Time History of Surge Velocity for Case 2. 
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SWAY VELOCITY 
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Figure 127. Time History of Sway Velocity for Case 2. 


Figures 128 and 129 show a plot of surge velocity (U), and sway velocity (V) 
versus time respectively for the x,,=-0.75 % L with a diveplane angle (és) of -20 
degrees, and rudder angle (é6r) of 0 degrees. Here we can see the periodic out of plane 
behavior of the steady state solution. This is what we expected from the eigenvalue 
analysis that we did in the previous section. 

Figures 130 and 131 show a plot of surge velocity (U), and sway velocity (V) 
versus time respectively for x,,=0 % L with a dive plane angle (és) of -20 degrees, and 
rudder angle (é6r) of 0 degrees. And once again, the steady state results of the numerical 


integration method match the linearized dynamic results exactly. 
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Figure 128. Time History of Surge Velocity for Case 3. 
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Figure 129. Time History of Sway Velocity for Case 3. 
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Figure 130. Time History of Surge Velocity for Case 4. 
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Figure 131. Time History of Sway Velocity for Case 4. 
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X. CONCLUSIONS AND RECOMMENDATIONS 


The problem of steady state response and dynamic stability analysis of submarines 


in free positive buoyancy ascent has been studied. The main conclusions of this work can 


be summarized as follows: 


Steady state motion is not, in general, restricted to the vertical plane, there can 
exist complicated out-of-plane motions. A combination of analytical results for the 
in plane motions and numerical simulations and continuation was used in order to 
identify all possible steady states. 


For some steady state solutions there can be pitchfork and hysteresis bifurcations. 
In fact, it was shown that the pitchfork is the primary mechanism for generating 
out of plane motions. 


The response of the vehicle can be very sensitive to system parameters. In such a 
case, any numerical integration results should be viewed with extreme caution. 


The bifurcation graphs that we generated are very useful in order to understand the 
steady state solution set of the vehicle for any parameter combination. 


As a recommendation for future research, continuation methods should be used to 
investigate the nature of the Hopf bifurcation, specifically stability of periodic 
solutions. 


Furthermore, the effect of nonzero propulsion should be investigated. It is possible 
that small nonzero propulsion forces could affect significantly the bifurcations. 
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APPENDIX 


1. BIFPACK 
Bifpack is a program package for calculating bifurcations wnitten in 
FORTRAN by R. Seydel [Ref. 15] 
In this work we used Bifpack to find the bifurcation characteristics of steady 
state solutions of submersibles. The resulting set of algebraic equations was solved by 


using PACKA which is devoted to system of nonlinear equations 
my, -) =0 


Here y and f(y,A) are vectors with n components. 

To use Bifpack we did not bother about how to program our original f(y,A) 
into the enlarged form that is used by Bifpack. This larger frame that embraces the 
original n scalar algebraic equations by Bifpack is called FRAMEA in Bifpack. To solve 
our problem we duplicated this frame and gave it different a name. This became our 
main program and then we entered n scalar formulas f,, f, ,f;,..... , f, into subroutine 
FCN. A sample of those two programs is included in the following sections. In 
conjunction with subroutine FCN, we also wrote functions which compute the integral 
values in Equations 8, 9, 11, and 12. These are functions FI2, FI3, FI5, FI6 in 


subroutine FCN. 


119 


OO@0O 


Oma@o 


@ O) @ 


2. MAIN PROGRAM 


SUB2.FOR 
PROGRAM MAIN 
IMPLICIT DOUBLE PRECISION (A-H,O-Z) 


DOUBLE PRECISION EPS 


EXTERNAL FCN 


INTERACTIVE WORK WITH THE PACK-A SUBPACKAGE OF 
BIFPACK 


COMMON/DATA/DS,DB,DELB,XGB,ZGB, XB,ZB,DRS,DRB 
OPEN (11, FILE = ’DATAAI1’,STATUS=’OLD’) 

OPEN (8 , FILE = DIAGRAM’ STATUS =’OLD’) 

OPEN (3 , FILE = ’START’ STATUS =’OLD’) 

OPEN (1 , FILE = ’OPTIONS’ SSTATUS=’OLD’) 

OPEN (10, FILE =’FORTRAN.MAT’ STATUS=’OLD’) 
OPEN (20, FILE =’FORTRANI.MAT’ STATUS =’OLD’) 
OPEN (15, FILE=’DEN.MAT’ STATUS =’OLD’) 

OPEN (16, FILE=’DEN1.MAT’ STATUS =’OLD’) 


PROBLEM DEPENDENT VARIABLES: 
EPS = (REL.) ACCURACY, N = NUMBER OF EQUATIONS 


EPS =1.E-5 
N=9 
JACOBI=1 


DATA 


PI =4.0*ATAN(1.0) 
WEIGHT = 12000.0 
L=17.425 
DB=0.0*PI/180.0 
DELB=2.0*WEIGHT/100.0 
XGB=-1.0*L/100.0 
ZGB=0.1 

XB=0.0*L/100.0 

ZB=0.0 

DRS =0.0*PI/180.0 
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DRB=0.0*PI/180.0 

WRITE (*,*) "WHAT IS DRS’ 
READ (*,*) DRS 

DRS =DRS*PI/180.0 

WRITE (*,*) DRS 


WRITE (11,*) DS AS PARAMETER ZGB=.1,DB IS 2 PERCENT 
OF WEIGHT’ 


INITIAL SOLUTION USED FOR THIS RUN 


Y(1)=8.6 U 
Y(2)=0 V 
Y(3)=-0.25 W 
Y(4)=0 Pp 
Y(5)=0 Q 
Y(6)=0 R 
Y(7)=0 PHI 
Y(8)=0.53 THETA 
Y(9)=0 PSIDOT 


WRITE (*,*) ’ Previous data in DIAGRAM to be purged?’ 
WRITE (*,*)’ ENTER O FOR PURGING, ENTER 1 FOR 


APPENDING:’ 


23 
29 


READ (*,*) NPURGE 
IF (NPURGE.EQ.0) REWIND 8 
IF (NPURGE.EQ.1) THEN 

DO 23 I=1,9999 

READ (8,’( )’, END=29) 

CONTINUE 

END IF 


CALL INTERA (N,FCN,EPS JACOBI) 

WRITE (*,*) ’ PRINT-OUTPUT IN FILE "DATAAI1"’ 
STOP 

END 


3. SUBROUTINE FCN. 


EOM1.FOR 
SUBROUTINE FCN (TDUMMY,Y,F) 
IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
DOUBLE PRECISION Y,F,PAR,A,ETACO,ZETACO,PARCO, TDUMMY 
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DOUBLE PRECISION DB, DR, XG, ZG, YG, XPROP, I2, I3, YB, 
KPROP, 15, XB, ZB, BOY, I6 ,NPROP, U, V, W, 
FI2, FI3,FI5,FI6 

DOUBLE PRECISION WEIGHT, IX, IY, IZ, IXY, IYZ, IXZ, L, RHO, 
G, AO, CDO, M, XPP, XQQ, XRR, XPR, XUDOT, 
XWQ, XVP, XVR, XQDS, XQDB, XDRDR, XRDR, XVV, 
XWW,XVDR,XWDS,XWDB,XDSDS,XDBDB,XRES, YPDOT, 
YRDOT,YPQ, YOR, YVDOT, YP, YR, YVQ,YWD, YWR, YV, 
YVW, YDR, YDRB, ZVV, ZDS, ZDB, KPDOT, KVW, 
MQDOT, MPP, MPR, ZQDOT, ZPP, ZPR, ZRR, 
ZWDOT,ZQ,ZVP,ZVR,ZW,NRDOT,NPQ,NOR, KRDOT, 
KPQ, KQR, KVDOT, KP, KR, KVQ, KWP, KWR, KV, 
MRR, MWDOT, MQ, MVP, MVR, MW, MVV, MDS, MDB, 
NPDOT,NVDOT, NP, NR, NVQ, NWP, NWR, NV, NVwW, 
NDR, NDRB, XRDRB, XRDRS 

PARAMETER (NDIM=12,NDIML=13) 

DIMENSION Y(NDIML),F(NDIML) 


FOLLOWING THREE STATEMENTS ARE FOR CONTINUATION AN 
BRANCH SWITCHING PURPOSES: 


COMMON/DATA/DS ,DB, DELB,XGB,ZGB,XB,ZB,DRS,DRB 
COMMON /CONTRO/ ETACO,ZETACO,PARCO,NCO,N1CO,N2CO, 
KCO,IND2CO,NFLAG 
PAR=Y(N1CO) 
F(N1CO) =Y(KCO)-ZETACO*Y(IND2CO)-ETACO 


GEOMETRIC PROPERTIES 


WEIGHT =12000.0 
IX =1760.0 
TY =9450.0 
IZ=10700.0 
IXY =0.0 
TYZ=0.0 
IXZ=0.0 
L=17.425 
RHO =1.94 
G=32.2 
YG=0.1 
YB=0.0 
A0=2.0 
CD0=0.0057 
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MASS =WEIGHT/G 
SURGE HYDRODYNAMIC COEFFICIENTS 


XPP = 7.030E-03*0.5*RHO*L**4 
XQQ =-1.470E-02*0.5*RHO*L**4 
XRR = 4.010E-03*0.5*RHO*L**4 
XPR = 7.640E-04*0.5*RHO*L**4 
XUDOT =-7.580E-03*0.5*RHO*L**3 
XWQ =-1.920E-01*0.5*RHO*L**3 
XVP =-3.240E-03*0.5*RHO*L**3 
XVR = 1.890E-02*0.5*RHO*L**3 
XQDS = 2.610B-02*0.5*RHO*L**3 
XQDB =-2.600E-03*0.5*RHO*L**3 
XRDR =-8.180E-04*0.5*RHO*L**3 
XVV_ = 5.290B-02*0.5*RHO*L**2 
XWW = 1.710B-01*0.5*RHO*L**2 
XVDR = 1.730E-03*0.5*RHO*L**2 
XWDS = 4.600E-02*0.5*RHO*L**2 
XWDB = 9.660E-03*0.5*RHO*L**2 
XDSDS =-1.160E-02*0.5*RHO*L**2 
XDBDB =-8.070E-03*0.5*RHO*L**2 
XDRDR =-1.010E-02*0.5*RHO*L**2 
XRES = CD0*0.5*RHO*L**2 
XPROP = XRES*ALPHA**2 
XRDRB=XRDR 

XRDRS=XRDR 

XVDRS=XVDR 

XDRB=XVDR 


SWAY HYDRODYNAMIC COEFFICIENTS 


YPDOT = 1.270E-04*0.5*RHO*L**4 
YRDOT = 1.240E-03*0.5*RHO*L**4 
YPQ) = 4.125E-03*0.5*RHO*L**4 
YOR =-6.510E-03*0.5*RHO*L**4 
YVDOT =-5.550E-02*0.5*RHO*L**3 
YP = 3.055E-03*0.5*RHO*L**3 
YR = 2.970E-02*0.5*RHO*L**3 
YVQ_ = 2.360E-02*0.5*RHO*L**3 
YWP = 2.350E-01*0.5*RHO*L**3 
YWR = =-1.880E-02*0.5*RHO*L**3 
YV =-9.310E-02*0.5*RHO*L**2 
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YVW = 6.840E-02*0.5*RHO*L**2 
YDRS =+2.270E-02*0.5*RHO*L**2 
YDRB =+2.270E-02*0.5*RHO*L**2 


HEAVE HYDRODYNAMIC COEFFICIENTS 


ZQDOT =-6.810E-03*0.5*RHO*L**4 
ZPP = 1.270B-04*0.5*RHO*L**4 
ZPR = 6.670E-03*0.5*RHO*L**4 
ZRR =-7.350E-03*0.5*RHO*L**4 
ZWDOT =-2.430E-01*0.5*RHO*L**3 
ZQ =-1.350E-01*0.5*RHO*L**3 
ZVP =-4.810E-02*0.5*RHO*L**3 
ZVR = 4.550E-02*0.5*RHO*L**3 
ZW =-3.020E-01*0.5*RHO*L**2 
ZVV_ =-6.840E-02*0.5*RHO*L**2 
ZDS =-2.270E-02*0.5*RHO*L**2 
ZDB =-2.270E-02*0.5*RHO*L**2 


ROLL HYDRODYNAMIC COEFFICIENTS 


KPDOT =-1.010E-03*0.5*RHO*L**5 
KRDOT =-3.370E-05*0.5*RHO*L**5 
KPQ =-6.930E-05*0.5*RHO*L**5 
KQR = 1.680E-02*0.5*RHO*L**5 
KVDOT = 1.270E-04*0.5*RHO*L**4 
KP =-1.100E-02*0.5*RHO*L**4 
KR =-8.410E-04*0.5*RHO*L**4 


KVQ =-5.115E-03*0.5*RHO*L**4 
KWP =-1.270E-04*0.5*RHO*L**4 
KWR_ = 1.390E-02*0.5*RHO*L**4 


KV = 3.055E-03*0.5*RHO*L**3 
KVW =-1.870E-01*0.5*RHO*L**3 


PITCH HYDRODYNAMIC COEFFICIENTS 


MQDOT =-1.680E-02*0.5*RHO*L**5 
MPP 5.260E-05*0.5*RHO*L**5 
MPR 5.040E-03*0.5*RHO*L**5 
MRR =-2.860E-03*0.5*RHO*L**5 
MWDOT =-6.810E-02*0.5*RHO*L**4 
MQ =-6.860E-02*0.5*RHO*L**4 
MVP = 1.180E-03*0.5*RHO*L**4 
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MVR = 1.730E-02*0.5*RHO*L**4 
MW = 9.860E-02*0.5*RHO*L**3 
MVV =-2.510E-02*0.5*RHO*L**3 
MDS =-1.113E-02*0.5*RHO*L**3 
MDB = 1.113E-02*0.5*RHO*L**3 


YAW HYDRODYNAMIC COEFFICIENTS 


NPDOT =-3.370E-05*0.5*RHO*L**5 
NRDOT =-3.400E-03*0.5*RHO*L**5 
NPQ =-2.110E-02*0.5*RHO*L**5 

NOR = 2.750E-03*0.5*RHO*L**5 
NVDOT = 1.240E-03*0.5*RHO*L**4 
NP =-8.405E-04*0.5*RHO*L**4 

NR  =-1.640E-02*0.5*RHO*L**4 

NVQ =-9.990E-03*0.5*RHO*L**4 
NWP =-1.750E-02*0.5*RHO*L**4 
NWR = 7.350E-03*0.5*RHO*L**4 
NV  =-7.420E-03*0.5*RHO*L**3 

NVW =-2.670E-02*0.5*RHO*L**3 
NDRS =-1.113E-02*0.5*RHO*L**3 
NDRB =+1.113E-02*0.5*RHO*L**3 


VARIABLE INITILIZATION 


RPM =0.0 

BOY = WEIGHT+DELB 
XG=XB+xXGB 

ZG =ZB+ZGB 

ALPHA =0.0 


TDUMMY IS DUMMY 

THIS ROUTINE EVALUATES THE PROBLEM DEFINING FUNCTION 
INPUT IS Y(1),...,Y(N), AND PAR=PARAMETER NOW, F(1),.. 

.F(N) OF THE FUNCTION OF RIGHT HAND SIDE OF F(Y,PAR) ARE 
ENTERED:(SIX DEGREE FREEDOM EQUATION OF MOTION OF C 
SUBMARINE) 


U=Y(1) 
V=Y(2) 
W=Y(3) 
P=Y(4) 

Q=Y(5) 
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R=Y(6) 
PHI=Y(7) 

THETA =Y(8) 
PSIDOT=Y(9) 
12=FI2(V,W,Q,R) 
13=FI3(V,W,Q,R) 
15 =FI5(V,W,Q,R) 
16=FI6(V,W,Q,R) 
SWAY =I2 
HEAVE=1I3 
PITCH =15 

YAW =I6 


PARAMETER IS DS FOR THIS RUN 
DS=PAR 
SURGE FORCE 


F(1) = MASS*V*R-MASS*W*Q+MASS*XG*Q**2+MASS*XG*R**2- 
MASS*YG*P*Q-MASS*ZG*P*R+XPP*P**2 + KQQ*Q**2 +XRR*R**2+ 
XPR*P*R+XWOQ*W*Q+XVP*V*P+XVR*V*R+U*Q* 

(XQDS*DS +XQDB*DB) + 
U*R*(XRDRS*DRS + XRDRB*DRB) + XVV*V**2+XWW*W**2+U*V* 
(XVDRS*DRS+XDRB*DRB) +U*W*(XWDS*DS + XWDB*DB) + 
(XDSDS*DS**2 +XDBDB*DB**2 +XDRDR*(DRS**2+DRB**2))*U**2- 
(WEIGHT-BOY) * DSIN(THETA) +XPROP*RPM*RPM-XRES*U*U 


R> Re Re Re Re RR 


SWAY FORCE 


F(2) = -MASS*U*R-MASS*XG*P*Q +MASS*YG*R**2-MASS*ZG*Q*R+ 
YPQ*P*Q+ YOR*Q*R+YP*U*P+ YR*U*R+YVQ*V*Q+YWP*W*P+ 
YWR*W*R+YV*U*V +YVW*V*W + YDRS*U**2*DRS + YDRB*U**2* 
DRB+ (WEIGHT-BOY)* DCOS(THETA)*DSIN(PHI) +MASS*W*P+ 
MASS*YG*P**2+SWAY 


R> Re RR 


HEAVE FORCE 


F(3) = MASS*U*Q-MASS*V*P-MASS*XG*P*R-MASS*YG*Q*R+ 
MASS*ZG*P**2 +MASS*ZG*Q**2+ZPP*P**2+ZPR*P*R+ 
ZRR*R**2+ZO0*U*Q+ZVP*V*P+ZVR*V*R+ZW*U*W + ZVV*V "824+ 
HEAVE+ U**2*(ZDS*DS + ZDB*DB) + (WEIGHT-BOY)* 
DCOS(THETA)* DCOS(PHI) 


R & & & 


126 


ROLL MOMENT 


G)i@) iC) 


F(4) = -IZ*Q*R+IY*Q*R-IXY*P*R+1YZ*Q**2-IYZ*R**2 +IXZ*P*Q+ 
MASS*YG*U*Q-MASS*YG*V*P-MASS*ZG*W*P+KPQ*P*Q+ 
KQR*Q*R+ KP*U*P+KR*U*R+KVQ*V*Q+KWP*W*P+KWR* 
W*R+KV*U*V+ KVW*V*W+ 
(YG*WEIGHT-YB*BOY)*DCOS(THETA)*DCOS(PHI)-(ZG*WEIGHT- 
ZB*BOY)*DCOS(THETA)*DSIN(PHI) +MASS*ZG*U*R 


R RR Rk 


PITCH MOMENT 


QAaAND 


F(5) = -IX*P*R+1Z*P*R+IXY*Q*R-IYZ*P*Q-IXZ*P**2 +IXZ*R**2- 
MASS*XG*U*Q+MASS*XG*V*P+MASS*ZG*V*R-MASS*ZG*W*Q+ 
MPP*P**2 + MPR*P*R+MRR*R**2 +MQ*U*Q+MVP*V*P+MVR*V*R 
+MW*U*W+4+MVV*V**2 +U**2*(MDS*DS + MDB*DB)-(XG*WEIGHT- 
XB*BOY)*DCOS(THETA)*DCOS(PHD- 
(ZG*WEIGHT-ZB*BOY)*DSIN(THETA)+PITCH 


CRE & 


YAW MOMENT 


AAD 


F(6) = -IY*P*Q+IX*P*Q+IXY*P**2-IXY *Q**2+IYZ*P*R-IXZ*Q*R- 
MASS*XG*U*R+MASS*XG*W*P-MASS*YG*V*R+MASS*YG*W*Q+ 
NPQ*P*Q+NQR*Q*R+NP*U*P+NR*U*R+NVQ*V*Q+NWP*W*P+ 
NWR*W*R+NV*U*V+ NVW*V*W-+NDRS*U**2*DRS +NDRB* 
U**2*DRB + (XG*WEIGHT-XB*BOY)*DCOS(THETA)*DSIN(PHD + 
(YG*WEIGHT-YB*BOY)*DSIN(THETA)+ YAW 


RB Re & & 


KINEMATIC RELATIONS 


Gare) 


F(7) = P+PSIDOT*DSIN(THETA) 
F(8) = Q-PSIDOT*DSIN(PHI)*DCOS(THETA) 
F(9) = R-PSIDOT*DCOS(PHI)*DCOS(THETA) 
RETURN 
END 


4. SAMPLE RUN 


MARS $ mun sub2 
WHAT IS DRS 
0. 
0.0000000000000000E + 00 
Previous data in DIAGRAM to be purged? 
ENTER 0Q FOR PURGING, ENTER 1 FOR APPENDING: 
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0 
THIS IS THE MENU YOU CAN CHOOSE FROM: 
QO : ENTERING STARTING DATA 
1 : CONTINUATION: NO RESTRICTIONS 
2 : CONTINUATION: SETTING OF OPTIONS, THEN RUN 
3 : CONTINUATION, FOLLOWING PREVIOUS OPTIONS 
-1 : ENDING THIS SESSION 
ENTER CHOICE of MODUS: 
0 
ENTER COMPONENT Y{( 1): 
4.2 
ENTER COMPONENT Y{( 2): 
0 
ENTER COMPONENT Y‘( 3): 
-.7 
ENTER COMPONENT Y‘( 4): 
0 
ENTER COMPONENT Y‘( 5): 
0 
ENTER COMPONENT Y{( 6): 
0 
ENTER COMPONENT Y{( 7): 
0 
ENTER COMPONENT Y{( 8): 
.1745 
ENTER COMPONENT Y{( 9): 
0 
ENTER PARAMETER: 
345 
THIS IS THE MENU YOU CAN CHOOSE FROM: 
0 : ENTERING STARTING DATA 
1 : CONTINUATION: NO RESTRICTIONS 
2 : CONTINUATION: SETTING OF OPTIONS, THEN RUN 
3 : CONTINUATION, FOLLOWING PREVIOUS OPTIONS 
ENDING THIS SESSION 
ENTER CHOICE of MODUS: 
2 
CURRENT STARTING DATA AT PARAMETER: 0.3450000000000000 
THESE ARE THE CURRENT OPTIONS: 
No: 1 , step length: -0.17406E-01 
No: 2 , index of fixed component: -1 
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No: 3 , level of continuation: 0 

No: 4 , maximum number of continuation steps : 200 

No: 5 , step bounds: Predictor:0.100E+00 , Param.:0.100E+00 , rel.change:0.10 
No: 6 , stability analysis (yves=1, no=0): 0 

No: 7, window : -0.35000E+00 .It. parameter .It. 0.35000E+00 

No: 8 , target point: no final value is set. 


ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 

ENTER -1 TO RETURN TO MENU: 

] 

ENTER INITIAL STEPSIZE (IN REAL): 

-.Q05 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 

ENTER -1 TO RETURN TO MENU: 
Zz 

CHOOSE INDEX FOR INITIAL STEP: 

ENTER 0 FOR CONTIN. WITH RESPECT TO PARAMETER, 
ENTER K FOR FIXING K-TH COMPONENT: 

0 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 

ENTER -1 TO RETURN TO MENU: 

3 

CHOOSE AMONG THREE LEVELS OF STEP CONTROL: 
ENTER 2. FOR FREELY VARIABLE STEPS, 

1 FOR VARIABLE STEP BUT FIXED INDEX, 
0 FOR BOTH STEPLENGTH AND INDEX FIXED: 

0 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 

ENTER -1 TO RETURN TO MENU: 

4 

Enter maximum NUMBER OF CONTINUATION STEPS: 

100 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 

ENTER -1 TO RETURN TO MENU: 

-] 

THIS IS THE MENU YOU CAN CHOOSE FROM: 

Q : ENTERING STARTING DATA 
1 : CONTINUATION: NO RESTRICTIONS 
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2 : CONTINUATION: SETTING OF OPTIONS, THEN RUN 
3 : CONTINUATION, FOLLOWING PREVIOUS OPTIONS 
-1 : ENDING THIS SESSION 
ENTER CHOICE of MODUS: 
2 
CURRENT STARTING DATA AT PARAMETER: 0.3450000000000000 
THESE ARE THE CURRENT OPTIONS: 


No: 1 , step length: -0.5Q000E-02 

No: 2 , index of fixed component: 0 

No: 3 , level of continuation: 0O 

No: 4 , maximum number of continuation steps : 100 

No: 5 , step bounds: Predictor:0.100E+00 , Param.:0.100E+00 , rel.change:0.10 
No: 6 , stability analysis (yes=1, no=0): 0O 

No: 7, window : -0.35000E+00 .It. parameter .It. 0.35000E+00 

No: 8 , target point: no final value is set. 


ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 

ENTER -1 TO RETURN TO MENU: 

0 
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